<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Flux Variability analysis (FVA)</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S3 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S4 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
.S5 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span style=' font-weight: bold;'>Flux Variability analysis (FVA)</span></h1><div  class = 'S1'><span style=' font-weight: bold;'>Authors :  Vanja Vlasov, Systems Biochemistry Group, LCSB, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Reviewers : Anne Richelle, Lewis Lab at University of California, San Diego.</span></div><div  class = 'S1'><span>Flux variability analysis (FVA) is a widely used computational tool for evaluating the minimum and maximum range of each reaction flux that can still satisfy the constraints using a double LP problem (i.e. a maximization and a subsequent minimization) for each reaction of interest [1]. </span></div><div  class = 'S1'><span></span></div><div  class = 'S2'><span texencoding="\begin{array}{lll}
\\{v}_{j,upper}/\ {v}_{j,lower}=max\limits _{v}/\min\limits _{v} &amp; v_{j} \\
\text{s.t.} &amp; Sv=0,\\
 &amp; l\leq v\leq u\\
\end{array}" style="vertical-align:-40px"><img src="" width="238.5" height="92" /></span></div><div  class = 'S1'><span>where </span><span texencoding="v \in R^{n}" style="vertical-align:-5px"><img src="" width="43" height="19" /></span><span> is the vector of specific reaction rates (metabolic fluxes) and  </span><span texencoding="{v}_{j,upper}" style="vertical-align:-6px"><img src="" width="41" height="20" /></span><span>  and  </span><span texencoding="{v}_{j,lower}" style="vertical-align:-6px"><img src="" width="40.5" height="20" /></span><span> are respectively the upper and lower values of each flux </span><span texencoding="{v}_{j}" style="vertical-align:-6px"><img src="" width="14.5" height="20" /></span><span> satisfying the system of linear equations.</span></div><div  class = 'S1'><span>Depending on the size of the model you are using for the analysis, use:</span></div><ul  class = 'S3'><li  class = 'S4'><span style=' font-family: monospace;'>fluxVariability()</span><span> function - for the low dimensional FVA;</span></li><li  class = 'S4'><span style=' font-family: monospace;'>fastFVA()</span><span> function - for the models with more than 1,000 reactions;</span></li><li  class = 'S4'><a href = "https://github.com/opencobra/COBRA.jl"><span>distributedFBA.jl</span></a><span> - for high dimensional FVA,</span><span style=' font-weight: bold;'> </span><span>models larger than 10,000 reactions </span><span texencoding="$^2" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>.</span></li></ul><h2  class = 'S5'><span>EQUIPMENT SETUP</span></h2><div  class = 'S1'><span>If necessary, initialize the cobra toolbox</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: pre"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0F5A0DB2" data-testid="output_0" data-width="428" data-height="927" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">      _____   _____   _____   _____     _____     |
     /  ___| /  _  \ |  _  \ |  _  \   / ___ \    |   COnstraint-Based Reconstruction and Analysis
     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - 2021
     | |     | | | | |  _  { |  _  /  |  ___  |   |
     | |___  | |_| | | |_| | | | \ \  | |   | |   |   Documentation:
     \_____| \_____/ |_____/ |_|  \_\ |_|   |_|   |   <a href="http://opencobra.github.io/cobratoolbox" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">http://opencobra.github.io/cobratoolbox</a>
                                                  | 

 &gt; Checking if git is installed ...  Done (version: 2.20.1).
 &gt; Checking if the repository is tracked using git ...  Done.
 &gt; Checking if curl is installed ...  Done.
 &gt; Checking if remote can be reached ...  Done.
 &gt; Initializing and updating submodules (this may take a while)... Done.
 &gt; Adding all the files of The COBRA Toolbox ...  Done.
 &gt; Define CB map output... set to svg.
 &gt; TranslateSBML is installed and working properly.
 &gt; Configuring solver environment variables ...
   - [-*--] ILOG_CPLEX_PATH: /usr/local/bin/cplex/ibm/ILOG/CPLEX_Studio1210/cplex/matlab/x86-64_linux
   - [-*--] GUROBI_PATH: /usr/local/bin/gurobi811/linux64/matlab
   - [-*--] TOMLAB_PATH: /usr/local/bin/tomlab/tomsym
   - [----] MOSEK_PATH: --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   Done.
 &gt; Checking available solvers and solver interfaces ...No valid license for the SNOPT solver


tomRun: ERROR! snopt solver returned empty Result

Run startup to set paths


 Done.
 &gt; Setting default solvers ... Done.
 &gt; Saving the MATLAB path ... Done.
   - The MATLAB path was saved as ~/pathdef.m.

 &gt; Summary of available solvers and solver interfaces

			Support 	   LP 	 MILP 	   QP 	 MIQP 	  NLP
	----------------------------------------------------------------------
	gurobi       	active        	    1 	    1 	    1 	    1 	    -
	ibm_cplex    	active        	    1 	    1 	    1 	    1 	    -
	tomlab_cplex 	active        	    1 	    1 	    1 	    1 	    -
	glpk         	active        	    1 	    1 	    - 	    - 	    -
	mosek        	active        	    0 	    - 	    0 	    - 	    -
	matlab       	active        	    1 	    - 	    - 	    - 	    1
	pdco         	active        	    1 	    - 	    1 	    - 	    -
	quadMinos    	active        	    1 	    - 	    - 	    - 	    -
	dqqMinos     	active        	    1 	    - 	    1 	    - 	    -
	cplex_direct 	active        	    1 	    1 	    1 	    - 	    -
	cplexlp      	active        	    1 	    - 	    - 	    - 	    -
	qpng         	passive       	    - 	    - 	    1 	    - 	    -
	tomlab_snopt 	passive       	    - 	    - 	    - 	    - 	    1
	lp_solve     	legacy        	    1 	    - 	    - 	    - 	    -
	----------------------------------------------------------------------
	Total        	-             	   11 	    5 	    7 	    3 	    2

 + Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.


 &gt; You can solve LP problems using: 'gurobi' - 'ibm_cplex' - 'glpk' - 'matlab' - 'pdco' - 'quadMinos' - 'cplexlp' 
 &gt; You can solve MILP problems using: 'gurobi' - 'ibm_cplex' - 'glpk' 
 &gt; You can solve QP problems using: 'gurobi' - 'ibm_cplex' - 'pdco' 
 &gt; You can solve MIQP problems using: 'gurobi' - 'ibm_cplex' 
 &gt; You can solve NLP problems using: 'matlab' 

&gt; Checking for available updates ... skipped</div></div></div></div></div><div  class = 'S8'><span>For solving line</span><span>ar programming problems in FBA and FVA analysis, certain solvers are required. The present tutoria</span><span>l can ru</span><span>n with glpk package, which does not require additional installation and configuration. Although, for the analysis of large models is recommended to use the GUROBI package.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: pre"><span >changeCobraSolver (</span><span style="color: rgb(170, 4, 249);">'ibm_cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'LP'</span><span >);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="9A9D0A9E" data-testid="output_1" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path.
 &gt; ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2019a on your operating system.</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >changeCobraSolver (</span><span style="color: rgb(170, 4, 249);">'ibm_cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'QP'</span><span >);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="F5A3B3F4" data-testid="output_2" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path.
 &gt; ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2019a on your operating system.</div></div></div></div></div><h2  class = 'S5'><span>PROCEDURE</span></h2><div  class = 'S1'><span>In thi</span><span>s tutorial, we will use the generic model of the human cellular metabolism, Recon2.0 </span><span texencoding="$^3" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>. Load the model</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">global </span><span >CBTDIR</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'Recon2.0model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >model = readCbModel(modelFileName);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="027BE71C" data-testid="output_3" data-width="428" data-height="18" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Each model.subSystems{x} is a character array, and this format is retained.</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'></div></div></div><div  class = 'S8'><span>The metabolites structures an</span><span>d reactions in Recon2.0 can be found in the Virtual Metabolic Human database (VMH, </span><a href = "http://vmh.life"><span>http://vmh.life</span></a><span>).</span></div><div  class = 'S1'><span>Constrain the model to limit the availability of carbon and oxygen energy sources. Find the uptake exchange reactions using </span><span style=' font-style: italic;'>findExcRxns</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[selExc, selUpt] = findExcRxns(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >uptakes = model.rxns(selUpt);</span></span></div></div></div><div  class = 'S8'><span>Select from the set of reactions defined in </span><span style=' font-style: italic;'>uptakes</span><span> those which contain a least one carbon in the metabolites involved in the reaction and set the flux values of these reactions to zero:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >subuptakeModel = extractSubNetwork(model, uptakes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >hiCarbonRxns = findCarbonRxns(subuptakeModel,1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >modelalter = changeRxnBounds(model, hiCarbonRxns, 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>Also close the other reaction related to the exchange of oxygen and energy sources:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >energySources = {</span><span style="color: rgb(170, 4, 249);">'EX_adp'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_amp(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_atp(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_co2(e)'</span><span >;</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_coa(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_fad(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_fe2(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_fe3(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_gdp(e)'</span><span >;</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_gmp(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_gtp(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_h(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_h2o(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_h2o2(e)'</span><span >;</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_nad(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_nadp(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_no(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_no2(e)'</span><span >; </span><span style="color: rgb(170, 4, 249);">'EX_o2s(e)'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >modelalter = changeRxnBounds (modelalter, energySources, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>For this tutorial, we will analyse the variability of several reactions of the human cellular metabolism in aerobic and anaerobic condition. For each simulation, the original model will be copied to a new model structure (e.g., </span><span style=' font-style: italic;'>modelfva1</span><span> for aerobic condition and </span><span style=' font-style: italic;'>modelfva2</span><span> for anaerobic condition). This preserves the constraints of the original model and allows to perform simulations with new constraints. Additionally, this method of renaming the model avoids confusion while performing multiple simulations at the same time.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelfva1 represents aerobic condition</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelfva1 = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelfva1 = changeRxnBounds(modelfva1, </span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelfva1 = changeRxnBounds(modelfva1, </span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >, -1000, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelfva2 represents anaerobic condition</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelfva2 = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelfva2 = changeRxnBounds(modelfva2, </span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >modelfva2 = changeRxnBounds(modelfva2, </span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >,  0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><h2  class = 'S5'><span>1) Standard FVA</span></h2><div  class = 'S1'><span>The full spectrum of flux variability analysis options can be accessed using the command:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [minFlux, maxFlux, Vmin, Vmax] = fluxVariability(model,...</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% optPercentage,osenseStr, rxnNameList, verbFlag, allowLoops, method);</span></span></div></div></div><div  class = 'S8'><span>The </span><span style=' font-family: monospace;'>optPercentage</span><span> parameter allows one to choose whether to consider solutions that give at least a certain percentage of the optimal solution (default - 100). Setting the parameters  </span><span style=' font-family: monospace;'>osenseStr = 'min'</span><span> or </span><span style=' font-family: monospace;'>osenseStr = 'max'</span><span> determines whether the flux balance analysis problem is first solved with minimization or maximisation (default - 'max'). The </span><span style=' font-family: monospace;'>rxnNameList</span><span> accepts a cell array list of reactions to selectively perform flux variability upon (default - all reactions of the model). This is useful for high-dimensional models as computation of flux variability for all reactions can be time consuming:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Selecting several reactions of the model that we want to analyse with FVA</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >rxnsList = {</span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >; </span><span style="color: rgb(170, 4, 249);">'ACOAHi'</span><span >; </span><span style="color: rgb(170, 4, 249);">'ALCD21_D'</span><span >; </span><span style="color: rgb(170, 4, 249);">'LALDO'</span><span >; </span><span style="color: rgb(170, 4, 249);">'ME2m'</span><span >;</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'AKGDm'</span><span >; </span><span style="color: rgb(170, 4, 249);">'PGI'</span><span >; </span><span style="color: rgb(170, 4, 249);">'PGM'</span><span >; </span><span style="color: rgb(170, 4, 249);">'r0062'</span><span >};</span></span></div></div></div><div  class = 'S8'><span>The </span><span style=' font-family: monospace;'>verbFlag</span><span> input determines the verbose output (default - false). </span><span style=' font-family: monospace;'>allowLoops </span><span>input determines whether loops are allowed in the solution (default - true). Note that </span><span style=' font-family: monospace;'>allowLoops==false</span><span> invokes a mixed integer linear programming implementation of thermodynamically constrained flux variability analysis for each minimization or maximisation of a reaction rate. The </span><span style=' font-family: monospace;'>method</span><span> parameter input determines whether are the output flux vectors also minimise the </span><span style=' font-family: monospace;'>0-norm</span><span>, </span><span style=' font-family: monospace;'>1-norm</span><span> or </span><span style=' font-family: monospace;'>2-norm</span><span> whilst maximising or minimising the flux through one reaction (default - 2-norm). </span></div><div  class = 'S1'><span>Run </span><span style=' font-family: monospace;'>fluxVariability()</span><span> on both models (</span><span style=' font-family: monospace;'>modelfva1</span><span>, </span><span style=' font-family: monospace;'>modelfva2</span><span>) to generate the minimum and maximum flux values of selected reactions (</span><span style=' font-style: italic;'>rxnsList</span><span>) in the model.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Run FVA analysis for the model with the constraints that simulates aerobic conditions:</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >[minFlux1, maxFlux1, Vmin1, Vmax1] = fluxVariability(modelfva1, 100, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >, rxnsList)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="76F592F3" data-testid="output_4" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">OPTIMAL</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsErrorElement scrollableOutput" uid="329B1E78" data-testid="output_5" data-width="428" data-height="170" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-errorType" style="white-space: normal; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">Error using solveCobraQP (line 1049)<br style="white-space: pre-wrap; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">[gurobi] Primal optimality condition in solveCobraQP not satisfied, residual = 0.00038861, while feasTol = 1e-06</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;"><br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">Error in fluxVariability&gt;getMinNorm (line 707)<br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">        solution = solveCobraQP(LPproblem, solverVarargin.QP{:});<br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;"><br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">Error in fluxVariability&gt;calcSolForEntry (line 685)<br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">        V = getMinNorm(LPproblem, LPsolution, numel(model.rxns), Flux, model, method, solverVarargin);<br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;"><br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">Error in fluxVariability (line 557)<br style="white-space: pre; font-style: normal; color: rgb(230, 0, 0); font-size: 12px;">        [maxFlux(i), V] = calcSolForEntry(model, rxnID ,LPproblem, method, allowLoopsI, minNorm, solverVarargin, heuristicSolutions{preCompMaxSols(i)}, -1);</div></div></div></div></div></div><div  class = 'S8'><span></span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Run FVA analysis for the model with the constraints that </span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% simulates anaerobic conditions:</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[minFlux2, maxFlux2, Vmin2, Vmax2] = fluxVariability(modelfva2, [], [], rxnsList) </span></span></div></div></div><div  class = 'S8'><span>The additional </span><span style=' font-family: monospace;'>n × k</span><span> output matrices </span><span style=' font-family: monospace;'>Vmin</span><span> and </span><span style=' font-family: monospace;'>Vmax</span><span> return the flux vector for each of the </span><span style=' font-family: monospace;'>k ≤ n</span><span> fluxes selected for flux variability.</span></div><div  class = 'S1'><span>You can further plot and compare the FVA results for the selected reaction from both models:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ymax1 = maxFlux1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ymin1 = minFlux1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ymax2 = maxFlux2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ymin2 = minFlux2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >maxf = table(ymax1, ymax2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >minf = table(ymin1, ymin2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >maxfxs = table2cell(maxf);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >minfxs = table2cell(minf);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >figure</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >plot1 = bar(cell2mat(maxfxs(1:end, :)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >plot2 = bar(cell2mat(minfxs(1:end, :)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">off</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >xticklabels({</span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ACOAHi'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ALCD21__D'</span><span >, </span><span style="color: rgb(170, 4, 249);">'LALDO'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >             </span><span style="color: rgb(170, 4, 249);">'ME2m'</span><span >, </span><span style="color: rgb(170, 4, 249);">'AKGDm'</span><span >, </span><span style="color: rgb(170, 4, 249);">'PGI'</span><span >, </span><span style="color: rgb(170, 4, 249);">'PGM'</span><span >, </span><span style="color: rgb(170, 4, 249);">'r0062'</span><span >})</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >set(gca, </span><span style="color: rgb(170, 4, 249);">'XTickLabelRotation'</span><span >, -80);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >yticks([-1000 -800 -600 -400 -200 0 200 400 600 800 1000])</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Reactions from the models'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Fluxes'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >legend({</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >}, </span><span style="color: rgb(170, 4, 249);">'Location'</span><span >, </span><span style="color: rgb(170, 4, 249);">'southwest'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >title(</span><span style="color: rgb(170, 4, 249);">'Variations in fluxes in the aerobic and anaerobic conditions'</span><span >)</span></span></div></div></div><h2  class = 'S5'><span>2) Fast FVA</span></h2><div  class = 'S1'><span>When the number of reaction on which you want to perform a flux variability is higher to 1000, we recommend using </span><span style=' font-family: monospace;'>fastFVA()</span><span> function instead of </span><span style=' font-family: monospace;'>fluxVariability(). </span><span>Note that for large models the memory requirements may become prohibitive.</span></div><div  class = 'S1'><span>The </span><span style=' font-family: monospace;'>fastFVA()</span><span> function only supports the </span><a href = "https://opencobra.github.io/cobratoolbox/docs/solvers.html "><span>CPLEX</span></a><span> solver.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S15'><span style="white-space: pre"><span >changeCobraSolver (</span><span style="color: rgb(170, 4, 249);">'ibm_cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'all'</span><span >, 1);</span></span></div></div></div><div  class = 'S1'><span>Run fastFVA analysis for the whole model (i.e. flux varaibility analysis performed on all reactions included in the model) with the constraints that simulates aerobic conditions:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S15'><span style="white-space: pre"><span >[minFluxF1, maxFluxF1, optsol, ret, fbasol, fvamin, fvamax, statussolmin, statussolmax] = fastFVA(modelfva1);</span></span></div></div></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S15'></div></div></div><div  class = 'S1'><span>Run fast FVA analysis for the whole model with the constraints that simulates anaerobic conditions:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[minFluxF2, maxFluxF2, optsol2, ret2, fbasol2, fvamin2, fvamax2,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    statussolmin2, statussolmax2] = fastFVA(modelfva2);</span></span></div></div></div><div  class = 'S1'><span>Plot the results of the fast FVA and compare them between the aerobic and anaerobic models:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ymaxf1 = maxFluxF1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >yminf1 = minFluxF1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ymaxf2 = maxFluxF2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >yminf2 = minFluxF2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >maxf =table(ymaxf1, ymaxf2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >minf =table(yminf1, yminf2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >maxf = table2cell(maxf);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >minf = table2cell(minf);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >figure</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >plot3 = bar(cell2mat(maxf(1:end, :)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >plot4 = bar(cell2mat(minf(1:end, :)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">off</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >xticks([0 2000 4000 6000 8000 10600])</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >yticks([-1000 -800 -600 -400 -200 0 200 400 600 800 1000])</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'All reactions in the model'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Fluxes'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >legend({</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >})</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >title(</span><span style="color: rgb(170, 4, 249);">'Variations in fluxes in the aerobic and anaerobic conditions'</span><span >)</span></span></div></div></div><h2  class = 'S5'><span>REFERENCES </span></h2><div  class = 'S1'><span>[1] Gudmundsson, S., Thiele, I. Computationally efficient flux variability analysis. </span><span style=' font-style: italic;'>BMC Bioinformatics. </span><span>11, 489 (2010).</span></div><div  class = 'S1'><span>[2] Heirendt, L., Thiele, I., Fleming, R.M. DistributedFBA.jl: high-level, high-performance flux balance analysis in Julia. </span><span style=' font-style: italic;'>Bioinformatics.</span><span> 33 (9), 1421-1423 (2017).</span></div><div  class = 'S1'><span>[3] Thiele, I., et al. </span><span>A community-driven global reconstruction of human metabolism. </span><span style=' font-style: italic;'>Nat. Biotechnol., </span><span>31(5), 419–425 (2013).</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% *Flux Variability analysis (FVA)*
% *Authors :  Vanja Vlasov, Systems Biochemistry Group, LCSB, University of 
% Luxembourg.*
% 
% *Reviewers : Anne Richelle, Lewis Lab at University of California, San Diego.*
% 
% Flux variability analysis (FVA) is a widely used computational tool for evaluating 
% the minimum and maximum range of each reaction flux that can still satisfy the 
% constraints using a double LP problem (i.e. a maximization and a subsequent 
% minimization) for each reaction of interest [1]. 
% 
% 
% 
% $$\begin{array}{lll}\\{v}_{j,upper}/\ {v}_{j,lower}=max\limits _{v}/\min\limits 
% _{v} & v_{j} \\\text{s.t.} & Sv=0,\\ & l\leq v\leq u\\\end{array}$$
% 
% where $v \in R^{n}$ is the vector of specific reaction rates (metabolic fluxes) 
% and  ${v}_{j,upper}$  and  ${v}_{j,lower}$ are respectively the upper and lower 
% values of each flux ${v}_{j}$ satisfying the system of linear equations.
% 
% Depending on the size of the model you are using for the analysis, use:
%% 
% * |fluxVariability()| function - for the low dimensional FVA;
% * |fastFVA()| function - for the models with more than 1,000 reactions;
% * <https://github.com/opencobra/COBRA.jl distributedFBA.jl> - for high dimensional 
% FVA, models larger than 10,000 reactions $$^2$.
%% EQUIPMENT SETUP
% If necessary, initialize the cobra toolbox

initCobraToolbox(false) % false, as we don't want to update
%% 
% For solving linear programming problems in FBA and FVA analysis, certain solvers 
% are required. The present tutorial can run with glpk package, which does not 
% require additional installation and configuration. Although, for the analysis 
% of large models is recommended to use the GUROBI package.

changeCobraSolver ('ibm_cplex', 'LP');
changeCobraSolver ('ibm_cplex', 'QP');
%% PROCEDURE
% In this tutorial, we will use the generic model of the human cellular metabolism, 
% Recon2.0 $$^3$. Load the model

global CBTDIR
modelFileName = 'Recon2.0model.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);

%% 
% The metabolites structures and reactions in Recon2.0 can be found in the Virtual 
% Metabolic Human database (VMH, <http://vmh.life http://vmh.life>).
%% 
% Constrain the model to limit the availability of carbon and oxygen energy 
% sources. Find the uptake exchange reactions using _findExcRxns_

[selExc, selUpt] = findExcRxns(model);
uptakes = model.rxns(selUpt);
%% 
% Select from the set of reactions defined in _uptakes_ those which contain 
% a least one carbon in the metabolites involved in the reaction and set the flux 
% values of these reactions to zero:

subuptakeModel = extractSubNetwork(model, uptakes);
hiCarbonRxns = findCarbonRxns(subuptakeModel,1);
modelalter = changeRxnBounds(model, hiCarbonRxns, 0, 'b');
%% 
% Also close the other reaction related to the exchange of oxygen and energy 
% sources:

energySources = {'EX_adp'; 'EX_amp(e)'; 'EX_atp(e)'; 'EX_co2(e)';...
    'EX_coa(e)'; 'EX_fad(e)'; 'EX_fe2(e)'; 'EX_fe3(e)'; 'EX_gdp(e)';...
    'EX_gmp(e)'; 'EX_gtp(e)'; 'EX_h(e)'; 'EX_h2o(e)'; 'EX_h2o2(e)';...
    'EX_nad(e)'; 'EX_nadp(e)'; 'EX_no(e)'; 'EX_no2(e)'; 'EX_o2s(e)'};
modelalter = changeRxnBounds (modelalter, energySources, 0, 'l');
%% 
% For this tutorial, we will analyse the variability of several reactions of 
% the human cellular metabolism in aerobic and anaerobic condition. For each simulation, 
% the original model will be copied to a new model structure (e.g., _modelfva1_ 
% for aerobic condition and _modelfva2_ for anaerobic condition). This preserves 
% the constraints of the original model and allows to perform simulations with 
% new constraints. Additionally, this method of renaming the model avoids confusion 
% while performing multiple simulations at the same time.

% modelfva1 represents aerobic condition
modelfva1 = modelalter;
modelfva1 = changeRxnBounds(modelfva1, 'EX_glc(e)', -20, 'l');
modelfva1 = changeRxnBounds(modelfva1, 'EX_o2(e)', -1000, 'l');
% modelfva2 represents anaerobic condition
modelfva2 = modelalter;
modelfva2 = changeRxnBounds(modelfva2, 'EX_glc(e)', -20, 'l');
modelfva2 = changeRxnBounds(modelfva2, 'EX_o2(e)',  0, 'l');
%% 1) Standard FVA
% The full spectrum of flux variability analysis options can be accessed using 
% the command:

% [minFlux, maxFlux, Vmin, Vmax] = fluxVariability(model,...
% optPercentage,osenseStr, rxnNameList, verbFlag, allowLoops, method);
%% 
% The |optPercentage| parameter allows one to choose whether to consider solutions 
% that give at least a certain percentage of the optimal solution (default - 100). 
% Setting the parameters  |osenseStr = 'min'| or |osenseStr = 'max'| determines 
% whether the flux balance analysis problem is first solved with minimization 
% or maximisation (default - 'max'). The |rxnNameList| accepts a cell array list 
% of reactions to selectively perform flux variability upon (default - all reactions 
% of the model). This is useful for high-dimensional models as computation of 
% flux variability for all reactions can be time consuming:

% Selecting several reactions of the model that we want to analyse with FVA
rxnsList = {'DM_atp_c_'; 'ACOAHi'; 'ALCD21_D'; 'LALDO'; 'ME2m';...
    'AKGDm'; 'PGI'; 'PGM'; 'r0062'};
%% 
% The |verbFlag| input determines the verbose output (default - false). |allowLoops| 
% input determines whether loops are allowed in the solution (default - true). 
% Note that |allowLoops==false| invokes a mixed integer linear programming implementation 
% of thermodynamically constrained flux variability analysis for each minimization 
% or maximisation of a reaction rate. The |method| parameter input determines 
% whether are the output flux vectors also minimise the |0-norm|, |1-norm| or 
% |2-norm| whilst maximising or minimising the flux through one reaction (default 
% - 2-norm). 
%% 
% Run |fluxVariability()| on both models (|modelfva1|, |modelfva2|) to generate 
% the minimum and maximum flux values of selected reactions (_rxnsList_) in the 
% model.

% Run FVA analysis for the model with the constraints that simulates aerobic conditions:
[minFlux1, maxFlux1, Vmin1, Vmax1] = fluxVariability(modelfva1, 100, 'max', rxnsList)
%% 
% 

% Run FVA analysis for the model with the constraints that 
% simulates anaerobic conditions:
[minFlux2, maxFlux2, Vmin2, Vmax2] = fluxVariability(modelfva2, [], [], rxnsList) 
%% 
% The additional |n × k| output matrices |Vmin| and |Vmax| return the flux vector 
% for each of the |k ≤ n| fluxes selected for flux variability.
%% 
% You can further plot and compare the FVA results for the selected reaction 
% from both models:

ymax1 = maxFlux1;
ymin1 = minFlux1;
ymax2 = maxFlux2;
ymin2 = minFlux2;

maxf = table(ymax1, ymax2)
minf = table(ymin1, ymin2)
maxfxs = table2cell(maxf);
minfxs = table2cell(minf);

figure
plot1 = bar(cell2mat(maxfxs(1:end, :)));
hold on
plot2 = bar(cell2mat(minfxs(1:end, :)));
hold off
xticklabels({'DM_atp_c_', 'ACOAHi', 'ALCD21__D', 'LALDO',...
             'ME2m', 'AKGDm', 'PGI', 'PGM', 'r0062'})
set(gca, 'XTickLabelRotation', -80);
yticks([-1000 -800 -600 -400 -200 0 200 400 600 800 1000])
xlabel('Reactions from the models')
ylabel('Fluxes')
legend({'Aerobic', 'Anaerobic'}, 'Location', 'southwest')
title('Variations in fluxes in the aerobic and anaerobic conditions')
%% 2) Fast FVA
% When the number of reaction on which you want to perform a flux variability 
% is higher to 1000, we recommend using |fastFVA()| function instead of |fluxVariability().| 
% Note that for large models the memory requirements may become prohibitive.
% 
% The |fastFVA()| function only supports the <https://opencobra.github.io/cobratoolbox/docs/solvers.html  
% CPLEX> solver.

changeCobraSolver ('ibm_cplex', 'all', 1);
%% 
% Run fastFVA analysis for the whole model (i.e. flux varaibility analysis performed 
% on all reactions included in the model) with the constraints that simulates 
% aerobic conditions:

[minFluxF1, maxFluxF1, optsol, ret, fbasol, fvamin, fvamax, statussolmin, statussolmax] = fastFVA(modelfva1);
%%

%% 
% Run fast FVA analysis for the whole model with the constraints that simulates 
% anaerobic conditions:

[minFluxF2, maxFluxF2, optsol2, ret2, fbasol2, fvamin2, fvamax2,...
    statussolmin2, statussolmax2] = fastFVA(modelfva2);
%% 
% Plot the results of the fast FVA and compare them between the aerobic and 
% anaerobic models:

ymaxf1 = maxFluxF1;
yminf1 = minFluxF1;
ymaxf2 = maxFluxF2;
yminf2 = minFluxF2;

maxf =table(ymaxf1, ymaxf2);
minf =table(yminf1, yminf2);

maxf = table2cell(maxf);
minf = table2cell(minf);

figure
plot3 = bar(cell2mat(maxf(1:end, :)));
hold on
plot4 = bar(cell2mat(minf(1:end, :)));
hold off
xticks([0 2000 4000 6000 8000 10600])
yticks([-1000 -800 -600 -400 -200 0 200 400 600 800 1000])
xlabel('All reactions in the model')
ylabel('Fluxes')
legend({'Aerobic', 'Anaerobic'})
title('Variations in fluxes in the aerobic and anaerobic conditions')
%% REFERENCES 
% [1] Gudmundsson, S., Thiele, I. Computationally efficient flux variability 
% analysis. _BMC Bioinformatics._ 11, 489 (2010).
% 
% [2] Heirendt, L., Thiele, I., Fleming, R.M. DistributedFBA.jl: high-level, 
% high-performance flux balance analysis in Julia. _Bioinformatics._ 33 (9), 1421-1423 
% (2017).
% 
% [3] Thiele, I., et al. A community-driven global reconstruction of human metabolism. 
% _Nat. Biotechnol.,_ 31(5), 419–425 (2013).
##### SOURCE END #####
-->
</div></body></html>